Density wave patterns for fermionic dipolar molecules on a square optical lattice: 

Mean-field-theory analysis 



K. Mikelsons and J. K. Freericks 
Department of Physics, Georgetown University, Washington, DC, 20057 USA 
(Dated: January 11, 2013) 

We model a system of ultracold fermionic dipolar molecules on a two-dimensional square lat- 
tice. Assuming that the molecules are in their nondegenerate hyperfine ground state, and that the 
dipole moment is polarized perpendicular to the plane (as in the recent experiments on 40 K- S7 Rb 
molecules), we approximate these molecules as spinless fermions with long range repulsive dipolar 
interactions. We use mean field theory to obtain the restricted phase diagram as a function of the 
filling, the strength of interaction and the temperature. We find a number of ordered density wave 
phases in the system, as well as phase separation between these phases. A Monte Carlo analysis 
shows that the higher-period phases are usually suppressed in the exact solution. 
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I. INTRODUCTION 

Recent experimental work [lj in ultracold molecule 
physics has formed fermionic dipolar molecules of 40 K 
and 87 Rb in their rovibrational and hyperfine ground 
state. Quantum degenerate dipolar systems have long 
been sought after because their long-range interactions 
makes them interesting for quantum computing applica- 
tions and for understanding ordered phases that can arise 
from long range effects. These systems are both compli- 
cated, but tunable, because the long range, anisotropic 
dipolar interaction can be manipulated with an external 
electric field [|| , and should give rise to rich physics with 
variety of different ordered phases of matter Q . 

The dipolar interaction is attractive if the dipole mo- 
ments are aligned head to tail, and this can lead to un- 
desirable recombination effects if the molecules have an 
exothermic reaction possible, as with 40 K- 87 Rb molecules 
This can be mitigated by confining molecules to a 
narrow two-dimensional layer, and setting a strong elec- 
tric field perpendicular to this layer, so that the dipole 
moments of molecules are aligned with the external field 
and always align head-to-head and tail-to-tail. In this 
case, the inter-molecule interactions are always repul- 
sive. Note that wc assume we can create a single isolated 
plane of dipoles, so there is no possibility to form super- 
fluid pairing between planes, where the interactions are 
attractive again 

We assume the dipolar molecules are loaded onto a 
single square optical lattice with the molecules fully po- 
larized perpendicular to the plane. Since the 40 K- 87 Rb 
molecules can be prepared in their lowest hyperfine nu- 
clear spin state, the internal spin degree of freedom is ef- 
fectively frozen out, and the molecules can be described 
as spinless fermions. 

Previous studies of spinless fcrmion models 
mostly focused on short range interactions and used 
the simplification of going to the limit of large spatial 
dimensions. These works show a rich phase diagram 
with charge-density-wave ordering and phase separation. 
There also has been some work in two dimensions focused 



on stripe physics @. 

More recent work, carried out in the context of ul- 
tracold dipolar molecules [H, fToMT2l | . focused primarily 
on finding new exotic phases (liquid crystal, smectic, or 
ncmatic), often achieved by tuning the direction of the 
external field with respect to the plane of molecules, but 
has not elucidated the details of the full phase diagram 
for the density-wave phases. 

The paper is organized as follows: First, in Sec. lU wc 
introduce the model and describe the possible density- 
wave orderings. Sec. IHII provides the details of the mean- 
field-theory formalism and the numerical calculations. 
The results are presented in Sec. lIVi and they are further 
discussed in Sec. El which arc followed by the conclusions 
in Sec. [Vll 



II. MODEL 

We describe the rovibrational ground-state dipolar 
molecules in their lowest hyperfine state with the follow- 
ing Hamiltonian: 



H = H k + H p = -tJ2(4 



c, + h.c. + 



«) 



where cj(ci) is the fermion creation (annihilation) oper- 
ator for a fermionic dipolar molecule at site i, rii = ejej 
is the number operator, t is the hopping amplitude be- 
tween the adjacent sites, and Uij = U/\fi — rj\ 3 is the 
long-range dipole interaction, with r, the position vector 
for the site i. The lattice spacing is taken to be equal to 
one. With this model, we have assumed that the hopping 
occurs only between nearest-neighbor sites, and that the 
molecules are always in the lowest level of the periodic 
optical lattice potential (this second assumption is not 
always true when the interaction between the molecules 
is too large). We have also assumed that the form for 
the dipole interaction is valid even when the distance be- 
tween molecules is only one lattice spacing. 
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Due to the bipartite nature of the square lattice, the 
Hamiltonian can be rewritten in a particle-hole sym- 
metric form and one can restrict to fillings that satisfy 
/ = (n) < 0.5. 

Since there is no spin degree of freedom in this model, 
it can only support molecule density wave spatial order. 
Such order is described by its unit cell — a parallelogram 
spanned by the two fundamental translation vectors for 
a given density wave order. The area of the unit cell 
(N c ) gives the number of sites with independent average 
molecule density which is modulated periodically as the 
unit cell is tiled across the lattice. The number of the in- 
dependent order parameters for this density wave is equal 
to N c — 1 (the average density, or filling, is not counted 
as an order parameter, and is fixed in the calculation). 
Since there is only a finite number of non-equivalent unit 
cells of a given area on a square lattice, there is only a 
finite number of possible density wave orders. We have 
considered all such orders that have up to four different 
order parameters (N c < 5) and hence this is a restricted 
phase diagram restricted to low-period phases only. The 
corresponding fundamental translation vectors and the 
unit cells for the twelve candidate phases are shown in 

Fig. m 
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FIG. 1: Left: Fundamental translation vectors (1x1,02) for 
possible density wave phases considered in this work, grouped 
by the area of the unit cell. Unit cell 1A describes a homoge- 
neous system without order, while 2B corresponds to checker- 
board order. Right: Unit cells for the density wave phases. 
Only density wave orders corresponding to unit cells with the 
solid outline were found to be stabilized in this study. 



III. METHODOLOGY 

We solve the model using mean field theory (MFT). 
This can be justified, since the the interaction is long 
range and consequently each site is effectively coupled 
to any other site, thereby increasing the effective dimen- 



sion of the mean field and decreasing the role of fluc- 
tuations (nevertheless, this is only an approximate so- 
lution). In fact, due to the absence of a local interac- 
tion, the MFT is equivalent to the dynamical mean-field 
theory (DMFT) approach, which becomes exact in the 
infinite-dimensional limit. The absence of a spin degree 
of freedom also implies that the model is in the Ising 
universality class, with a finite transition temperature in 
2D. Thus, while MFT is expected to overestimate the 
transition temperature, it is also expected to give the 
qualitatively correct phase diagram. Within MFT, the 
interaction part of Hamiltonian is approximated as: 



i(rij) + {rii}nj - (ni)(nj) , 



(2) 



which means that the correlations of the density fluctu- 
ations are neglected. This corresponds to the first order 
(Hartrce-Fock) self-consistent perturbation theory result 
and is expected to be accurate for small U /t. 

In the MFT approximation, the order parameter (n,) 
is a fixed parameter in Hamiltonian, and acts as a site 
dependent potential. The resultant MFT Hamiltonian 
is quadratic in the (c, c') operators and can be easily 
diagonalized. We generically work on a lattice with a 
large, but finite, number of lattice sites N and periodic 
boundary conditions. Further simplification comes from 
exploiting translational invariancc, which means that the 
Hamiltonian is block-diagonal in reciprocal space. The 
size of these blocks is N c , and it grows with the num- 
ber of the order parameters while the number of blocks, 
N/N c , grows with the number of discretization points in 
reciprocal space, N, which is the same as the number of 
real-space lattice sites included in the calculation. The 
Hamiltonian in reciprocal space becomes: 



H 



N/N c 

EE' 



Ki + k 



(e(fc + Ki) - n)5ij + Vi- 



-Kj +k ' 



(3) 

where e(k) = —2t(cosk x + cosk y ) is the tight-binding 
dispersion, ji is the chemical potential, Ki are the recip- 
rocal lattice points corresponding to the real-space ba- 
sis vectors of a given order (see Fig. [1]), and k denotes 
the discrete summation points in the corresponding first 
Brillouin zone (FBZ) (see Fig. [2]). The interaction part 
Vij = U(Ki - Kj){m{Ki - Kj)) where m(K z ) is the 
Fourier transform of the real space density (order pa- 



rameter): m(Ki) = Y] 



e iKir *m 



U(K l ) = U 

all 3^0 



and U(Ki) is: 



(4) 



Since the the dipole interaction decays sufficiently fast 
with increasing distance, the sum can be carried out for 
an infinite range of distance, giving more precise values 
for the effective interaction in /c-space. As seen from 
Eq. [3J k is a good quantum number, while the spatially 
modulated mean field in the interaction term Vu causes 
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molecule scattering among the N c X-points. The calcu- 
lation starts with a random guess for the order parameter 
(rii) (i = 1, . . . , N c — 1), which is put into the Hamiltonian 
(Eq.[3]), and the full set of eigenvalues (e^ Q ) and the cor- 
responding eigenvectors [i/>k a (Ki)] are obtained by solv- 
ing the eigenvalue equation (for all k and a = 1 . . . N c ): 

£ [(e(fc + Ki) + V i3 ] 1>- ka {Ki) = e^MKi) . 

i=i 

(5) 

For simplicity, we denote e< = tj. a and ipy = i>% a {Kj). 
The energy per site and the filling fraction are next cal- 
culated via 



1 N 1 N 



(6) 



i=l 



respectively, with n £i = l/(l + e^ e<_ ^' T ) the Fermi-Dirac 
distribution and T being the temperature. The mean- 
field entropy per site is 

1 N 

S = -^^[n ei m(n ei ) + (l-n ei )ln(l-n ei )] , (7) 

i=l 

and the free energy satisfies F 1 = E—TS. The momentum 
dependent density is: 



N e 

E 1 



(8) 



Calculating the expectation value of the average real 
space density then yields the new values of the order 
parameter: 



(m{K 3 )) 



1 



N 

'n E 

1=1 1 



E 



(9) 

which can be used to calculate the structure factor 
S(Q) = \m(Q)\ 2 . These updated order parameter val- 
ues are inserted back into the Hamiltonian in Eq. ([3]), 
and the iterative process is repeated until the order pa- 
rameter has converged to a fixed point. The chemical 
potential needs to be adjusted after every MFT itera- 
tion, so as to maintain constant filling. 

Solving the MFT equations represents finding a global 
minimum for the free energy functional in a multidimen- 
sional space. Since there can be multiple distinct stable 
solutions (corresponding to local minima in the free en- 
ergy), the calculation needs to be repeated for several 
random initial choices of the order parameter, and only 
the solution corresponding to global free energy minimum 
should be taken. More starting points are necessary for 
higher order phases, since the parameter space of solu- 
tions is larger. For phases with only one order parameter 
(2B), the MFT solution can be found by a successive bi- 
section method. Overall, we found that with 20 to 50 



initial random order parameters, in more than half of 
the cases, the calculation converged to the lowest free 
energy minimum. The number of required integration 
points N depends on the value of U/t: for U/t — > oo 
taking N = N c is sufficient, while for small U/t finer 
discretization is needed to account for the fc-dependcncc 
of dispersion e(k). For a given U/t, we started with a 
smaller number of integration points N , and increased it 
by a factor of four, until further increase did not produce 
a significant change in the phase diagram. 

The search for the MFT solutions is repeated for all 
possible density- wave orders (as described above), and 
only the solution with the lowest free energy is ultimately 
taken. 




FIG. 2: (color online) Reciprocal space for the 3B order. The 
FBZ is shown in blue in the center of square. The black 
circles denote K points (N c = 3), associated with the density 
wave order, while the gray circles are the discretization k- 
points used in the calculation. The open circles denote other 
possible values of K + k within reciprocal space. The dotted 
arrows denote other reciprocal space vectors for 3B order used 
to construct the FBZ. 



A simplification occurs in the limit U/t — > oo. In 
this case, the kinetic energy term in the Hamiltonian in 
Eq. (fTJ is negligible, and the MFT Hamiltonian becomes: 



H 



Uij [ni(rij) + (ni)nj - (10) 



with 



W^ii + const , 



W i = 2Y t U ij {n i ) . 



(11) 



(12) 
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The solution yields readily: 

(n l ) = (l + e^-^ T y 1 . (13) 

Equations (| 12[) and (|T3| need to be solved self- 
consistent ly to yield the MFT result. 

To account for the possibility of phase separation, we 
need to calculate the free energy for fixed U jt and T/t for 
full range of filling / and perform a Maxwell construction. 
A clear indicator of phase separation is that the chemical 
potential dependence on the filling becomes multivalued 
(since ^ = dF/df), as will be shown in the following 
Section. 
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IV. RESULTS 

We present results for the phase diagrams, obtained 
both with and without phase separation considered. To 
illustrate the procedure for detecting phase separation, 
we show results for the chemical potential at T = and 
U/t = 4 in Fig. H The MFT solution is particle-hole 
symmetric, so that n(f) + fx{l — f) = 2/i(0.5). How- 
ever, it shows unphysical behavior at some regions where 
the chemical potential decreases with increasing filling. 
The physical solution avoids such regions by allowing for 
phase separation, which consists of mixtures of two dif- 
ferent phases with the same chemical potential, at the 
same temperature and interaction U/t. As seen in Fig. [31 
phase separation occurs at fillings slightly above or below 
the commensurate values, and is accompanied by a gap 
in the density of states at the commensurate values of 
filling. 

The phase diagram for U/t = oo is shown in Fig. 2] 
In this case, the kinetic-energy part of the Hamiltonian 
is suppressed (t = 0), and the order at zero temperature 
extends to all values of the filling. The dominant order 
is the checkerboard pattern (2B), commensurate at half 
filling. At low temperature, higher-period orders emerge 
around the corresponding commensurate fillings. Yet 
low-period orders are more stable against thermal fluc- 
tuations, as indicated by their persistence up to higher 
temperatures. The transition between the high temper- 
ature homogeneous (1A) to checkerboard phase (2B) is 
continuous, as is the transition between the 2B and 4D 
phases. The remaining phase transitions, which occur 
between higher-period phases at lower temperatures, are 
discontinuous and can involve phase separation. 

For finite U = 4£ (see Fig. [5J , the zero-temperature or- 
der is suppressed for small values of filling, yet the phase 
diagram around half-filling is relatively unchanged, as 
compared to U/t = oo. Again, the checkerboard phase is 
dominant at higher temperatures. Higher-period phases 
(3B, 4C and 5B), although still confined to low temper- 
atures, appear to be more protracted in the range of fill- 
ing. The phase separation between these phases is also 
more prominent, but not qualitatively different than in 
the U/t = co case. An interesting feature is the reentrant 



FIG. 3: (color online) Chemical potential for U/t — 4 and 
T = 0. The red dashed line shows the MFT solution [which 
obeys particle-hole symmetry, fi(f) +/i(l — /) = 2/i(0.5)], but 
includes regions of unphysical behavior with dfi/df < 0. The 
blue line shows the physical solution, which includes regions 
of phase separation with d/i/df = 0, accompanied by corre- 
sponding gaps at commensurate fillings. Two shaded regions 
near / = 0.35 have the same area, according to Maxwell con- 
struction. 



behavior for this phase around / = 0.2 (evidenced by the 
backward curving of the stability region of the checker- 
board phase). This is presumably induced by the van 
Hove singularity, which can stabilize the checkerboard 
order with the help of thermal activation. 

The phase diagram for T = is shown in Fig. [5] De- 
creasing the interaction from U/t = oo causes ordering 
to disappear for low filling, as the kinetic energy effect 
is dominant in a dilute system. As the interaction is de- 
creased, the phases with filling near commensurate values 
are cut off successively. Before its disappearance, phase 
5B is significantly expanded in filling around U/t = 2.5, 
perhaps due to a better nesting of Fermi surface, relative 
to other phases. Below, but near half filling, decreasing 
interaction suppresses the 4D phase in favor of weak order 
in the 4B phase and eventually leaving only the checker- 
board order (2B). The phase separation is not confined to 
the immediate vicinity of the commensurate fillings, (also 
seen in Figs. @] and 0, as it can replace other phases for 
a range of fillings. 

The checkerboard ordering is the most favorable phase 
at half filling for any value of U /t and it persists down to 
U/t = 0, due to perfect nesting of the Fermi surface (see 
Fig. [7]). For large U/t, the transition temperature T c is 
linear in U , the largest energy scale. (Since there is only 
one spin species, there is no energy scale associated with 
spin exchange, such as J = 4t 2 /U in the Hubbard model.) 
For small U /t, T c shows exponentially activated behavior 
and can be fit with a semi-analytical form, obtained by 
approximating the non-interacting density of states with 
a logarithmic divergence representing the van Hove sin- 
gularity. Similarly, we find an exponentially contracted 
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FIG. 4: (color online) U = oo phase diagram without (top) 
and with (bottom) phase separation. Phase separation is 
shown as striped regions. The transitions between 4D and 
2B (as well as between 2B and 1A) phases are continuous. 
The rest of the transitions are discontinuous because of phase 
separation. 



FIG. 5: (color online) U = 4t phase diagram without (top) 
and with (bottom) phase separation. There is no order at 
low filling, yet for / > 0.25, the phase diagram qualitatively 
resembles the U/t — oo result. Note the reentrant behavior 
for the checkerboard phase (2B) around / = 0.2 and the phase 
separation between the 5C and homogeneous phases. 



range of filling near U = (see Fig. [6} for this phase. 
Thus, checkerboard ordering near U = 0, while possible 
due to perfect nesting, is strongly suppressed due to the 
proximity to the van Hove singularity. The gap in the 
density of states as well as the order parameter at T = 
show very similar dependence on U/t. 

We also report the entropy per particle, since this 
parameter is important in experiments. Calculated at 
T c for the checkerboard phase at / = 0.5, the entropy 
per particle shows similar behavior with U/t to that 
of T c . However, for U/t = oo, the ordering starts at 
S(T = T c ) = 2 In 2, clearly an overestimate due to the 
neglected density fluctuations in the MFT. For smaller 
U /t, a lower entropy state has to be reached to observe 
the checkerboard order. 

The most common experimental test for spatial order 
is Bragg diffraction of light by the molecules, which di- 
rectly yields the structure factor in reciprocal space and 
can be used to identify the order [l3j]. Such an experi- 
ment is the analog of an x-ray diffraction experiment on 



a condensed-matter system used to determine its crys- 
tal structure. Table U shows a comparison of different 
density-wave phases found in this study by listing the 
real space density patterns, reciprocal space points and 
structure factors for U/t = 4. Each of the phases have 
sufficiently different signatures, both in the location and 
the relative strength of the peaks of S(Q), making these 
phases easily identifiable through Bragg diffraction mea- 
surements. 

Another, less direct way to observe spatial order is 
time-of-flight measurements that use the expansion of the 
molecular cloud to infer the initial momentum distribu- 
tion of the molecules. However, the momentum distribu- 
tion data do not show prominent differences amongst the 
different ordered phases, so using this approach to detect 
the ordering is difficult. 
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FIG. 6: (color online) T = phase diagram without (top) and 
with (bottom) phase separation. Decreasing the interaction 
contracts the range of filling of the ordered phases and pro- 
gressively eliminates phases commensurate with low values of 
filling. Only the checkerboard phase survives down to U = 0. 
Phase separation replaces the 4D phase near / = 0.28 and 
/ = 0.36 for larger U /t. In parts of the phase diagram, 4C 
and 5C phases show phase separation with the homogeneous 
state. 



V. DISCUSSION 

The MFT results show a diverse phase diagram with 
a number of distinct density-wave phases. However, the 
MFT is an approximate method, which neglects fluctu- 
ations and favors ordered states. This brings to mind 
the question: how many of these density-wave phases are 
actually present in the exact solution? 

In the limit U/t = oo, the Hamiltonian essentially de- 
scribes a lattice gas model and can be solved numeri- 
cally without approximation using classical Monte Carlo 
methods. Preliminary results support the existence of the 
checkerboard (2B) and 4D phases, although strongly sup- 
pressed in temperature and filling, compared to the MFT 
result. Since the MFT neglects fluctuations, it is not sur- 
prising that the range of ordered phases is contracted in 
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FIG. 7: (color online) / = 0.5 phase diagram. For large 
U/t, the transition temperature satisfies T c = aU/2, where 
a = U[K = (7T,7r)] = 1.323. For small U/t, T c fits perfectly 
to a semianalytical form. The gap in the density of states (not 
shown) at / = 0.5 and T = shows similar dependence on 
U/t. The entropy per particle (dotted line) at T = T c reaches 
S = 21n2 at U/t = oo. 



a more precise calculation. However, some higher-period 
phases become more prominent in the MFT solution for 
intermediate values of U /t, so their presence cannot be 
entirely ruled out by this U/t = oo study. A more com- 
plete study will be the topic of further research. 

Also, the phase separation needs to be reviewed us- 
ing better approximation techniques. The presence of 
phase separation is tied to the existence of higher-period 
phases and the fact that the MFT strongly stabilizes 
these phases at the corresponding commensurate fillings. 
Thus, phase separation is promoted by the MFT just 
like the corresponding ordered phases are. It is possible 
that including even higher order phases (beyond 4 order 
parameters) in a MFT study would replace the phase- 
separated regions with higher-period ordered phases that 
are commensurate at the corresponding values of filling, 
yielding a devil's staircase-like structure. Another pos- 
sible alternative to phase separation is the formation of 
a Wigner crystal like state. The reentrant behavior of 
the checkerboard phase at U = At in the MFT results 
is very likely an artifact of this method and should be 
reexamined using more precise techniques. 

It is notable that the MFT results do not favor stripes 
that arc aligned with the lattice (corresponding to phases 
2A, 3A, 4A or 5A). We found that these phases are ei- 
ther not stable or they are replaced by other phases with 
lower free energy. For small U/t, good nesting of the 
Fermi surface appears to be important to lower the en- 
ergy of the system, while for large U /t, strong repulsive 
interactions between adjacent sites prevents formation of 
such striped phases. We also see no evidence for any 
nematic or smectic "liquid-crystal" phases. 

Another important question is, which of these density- 
wave phases can be reached and detected in experiment? 
The MFT results show large entropy per particle when 
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phase fc-space points structure factors / T/t 



TABLE I: Structure factors for different phases at represen- 
tative values of filling and temperature for U/t = 4. The first 
column shows the real space patterns with darker squares de- 
noting larger average density per site. The second column 
shows the reciprocal space points corresponding to a given 
density-wave order, followed by the corresponding structure 
factors at these points, S(K), relative to So = S(K = 0). 
Since S(—K) = S(K), only the structure factors for K x < 
are shown. The last two columns give the corresponding fill- 
ing and temperature, respectively. 

the checkerboard state develops at half filling for large in- 
teraction, but unfortunately the entropy is severely over- 



estimated in this study, as shown by preliminary Monte- 
Carlo results. The most easily observable features in ex- 
periment are the structure factors and the gap in the 
density of states at the commensurate values of filling. 
These features are most pronounced for the checkerboard 
phase. In a trap environment, a large gap at half filling 
should translate into large spatial regions of checkerboard 
order if the local density approximation is accurate. Of 
course, it is not clear if phases that are seen in a homo- 
geneous system are stabilized or destabilized by the trap 
environment. This is another question to be addressed 
in a future study. 



VI. CONCLUSIONS 

We have studied the ordering of ultracold diatomic 
gases on a two-dimensional square optical lattice. We 
have formulated a model, based on spinless fermions with 
long range dipole interactions, and solved it using mean 
field theory. We have found a number of possible density- 
wave orders that occur mostly around the corresponding 
commensurate values of filling and create gaps in the den- 
sity of states at these values of the filling. The phase 
separation is also often seen between different ordered 
phases in the MFT results. Overall, the checkerboard 
phase is most stable near half filling and has the high- 
est transition temperature. Due to the lack of an en- 
ergy scale associated with spin exchange, the transition 
temperatures increase with interaction strength. Differ- 
ent density wave phases show characteristic and distinct 
structure factors. We have discussed the validity of these 
results and the limitations of the MFT method used in 
this study. 
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